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We describe a flexible and modular delayed-feedback nonlinear oscillator that 
is capable of generating a wide range of dynamical behaviours, from periodic oscil- 
lations to high-dimensional chaos. The oscillator uses electrooptic modulation and 
fibre-optic transmission, with feedback and filtering implemented through real-time 
digital-signal processing. We consider two such oscillators that are coupled to one 
another, and we identify the conditions under which they will synchronize. By ex- 
amining the rates of divergence or convergence between two coupled oscillators, we 
quantify the maximum Lyapunov exponents or transverse Lyapunov exponents of 
the system, and we present an experimental method to determine these rates that 
does not require a mathematical model of the system. Finally, we demonstrate a 
new adaptive control method that keeps two oscillators synchronized even when 
the coupling between them is changing unpredictably. 
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1. Introduction 

Private communications, fast physical random number generators and spatiotem- 
porally distributed sensor networks have provided the context for possible new ap- 
plications of chaotic dynamical systems. A key requirement for such applications is 
the development of reliable and robust generators of chaotic waveforms with broad 
spectral bandwidths. By reliable, we mean that given parameters of the system, 
the dynamical properties are reproducible, both experimentally and theoretically. 
Robust implies that the system exhibits chaotic behaviour over a region of parame- 
ter space (i.e., with few periodic windows.) There have been several realizations of 
delayed-feedback optoelectronic oscillators that meet these criteria. Systems that 
can be configured for integrated optoelectronic fabrication and can function at fre- 
quency ranges from tens of GHz down to kHz may find applications in acoustic, 
biological, chemical, electromagnetic, and mechanical scenarios on nano- to macro- 
scopic spatial scales (Argyris et al. 2005, Kouomou et al. 2005, Uchida et al. 2005, 
Uchida et al. 2008, lUing et al. 2007, Reidler et al. 2009, Sorrentino & Ott 2008, 
2009, Cohen et al. 2008, Yousefi et al. 2008, Argyris et al. 2008). 
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Our focus in this paper is on the dynamics of delayed-feedback nonhnear os- 
cillators constructed from modular optoelectronic components. Delayed feedback 
enables such systems to generate a wide variety of waveforms, with differing degrees 
of complexity that depend on the parameters used. In particular, the time-delay, 
feedback strength, and filter parameters can be tuned to produce highly stable 
periodic waveforms (Yao & Malcki 1996) as well as complex waveforms that are 
characteristic of robust, high-dimensional chaos (Peil et al. 2009). 

In §2 we introduce the basic optoelectronic system and the delay-differential 
equations used for a continuous-time description of the dynamics (Kouomou et al. 
2005, Cohen et al. 2008). The intrinsic nonlinearity of the system arises from the 
integrated optical Mach-Zehnder modulator which changes the intensity of light 
transmitted depending on the cosine squared of a modulation voltage applied to its 
electrodes. We then chart the dynamical behaviour of the system, using bifurcation 
diagrams, as the feedback strength and delay time of the feedback loop are varied. 
The complexity of the waveforms generated is assessed by the Lyapunov dimension, 
and we illustrate the wide range of dynamics accessible. 

In §3 we first motivate and then show how to incorporate digital signal process- 
ing (DSP) capabilities in the delayed-feedback system. This enables precise, real 
time control of system parameters, such as the time-delay and filter characteristics, 
in a flexible manner well-suited for applications in communications and sensor net- 
works. The transition from continuous-time to discrete-time equations is outlined; 
the system is now governed by finite difference equations that describe its time 
evolution in terms of the system state as sampled at discrete times by an analogue 
to digital converter (Toomey et al. 2009). Even though our DSP implementation 
was aimed at kHz frequencies, such systems can be extended easily into the GHz 
range. 

The question of isochronal synchronization of these nonlinear oscillators (Fischer 
et al. 2006, Klein et al. 2006, Rogers-Dakin et al. 2006, Schwartz & Shaw 2007, 
Zhou & Roy 2007, Pranz et al. 2008) is central to possible applications in sensor 
networks (Sorrentino & Ott 2008, 2009). We thus consider coupled oscillators next 
in §4, where the many different configurations in which even two oscillators may 
be coupled are outlined. We then restrict ourselves to the schemes that we have 
explored in some detail. A diffusive-coupling scheme that allows the coupled systems 
to synchronize and retain the dynamical behaviour of the uncoupled systems is of 
particular interest. Several results on the dependence of synchronization error on 
coupling strength that have been obtained mathematically are verified through 
numerical simulations and tested experimentally. In particular, we emphasize that, 
in the experiments, noise and differences in nominally matched system parameters 
are unavoidable. We idenfity parameter regimes for the coupling strength where 
stable synchronization is observed. 

While the steady-state synchronization error is an important quantity to mea- 
sure with regards to sensor and communications applications, the transients to- 
wards synchrony and away from synchrony are important as well, and we study 
these in §5. The time scales for these transients set the limits on communication 
rates and detection of environmental perturbations. One may determine the max- 
imum Lyapunov exponent for a dynamical system by measurement of transients 
away from synchrony (Cohen et al. 2008). When a mathematical model is avail- 
able, one may predict the dynamics of an experimental system for several delay 
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times by performing data assimilation of experimental data using synchronization 
of the mathematical model to the data (Cohen et al. 2008, Marino et al. 2009, 
Quinn et al. 2009, So et al. 1994, Sorrentino & Ott, 2009b). Further, it is possi- 
ble to estimate distributions for finite-time Lyapunov exponents of the system. It 
should be noted that given two replicas of a dynamical system, one may estimate 
Lyapunov exponents from transients even when one does not have a mathematical 
model of the system. 

For applications of synchronized chaotic systems to sensor networks a novel 
adaptive synchronization approach has been recently conceived by Sorrentino & 
Ott (2008, 2009). When the coupHng channels between diffusively coupled chaotic 
dynamical systems serving as nodes of the network are perturbed at time scales 
slow compared to those of the chaotic fluctuations, they showed that it is possible 
to not only maintain synchrony between the systems. In the process of doing so, it 
is also possible to estimate and track the time- varying perturbations of the coupling 
strengths. In the illustrative case of two coupled systems, we have recently shown 
(Ravoori et al. 2009) that this scheme can be implemented experimentally. In §6 
we describe the scheme as implemented in the DSP based system described in §3, 
and we examine its effectiveness in maintaining synchrony and tracking the time- 
dependent perturbations of the coupling channel. 

In §7 we summarize our results and discuss future directions of research. 



2. Chaotic optoelectronic oscillator 

Figure 1(a) shows a diagram of the chaotic optoelectronic oscillator considered 
here, composed of a laser, electrooptic intensity modulator, photoreceiver and elec- 
trical filter, all connected together in a time-delayed feedback loop. This system 
was originally considered by Neyer and Voges (1986), who recognized its potential 
for bistability and chaos. The system was later adapted for use as a high-quality 
microwave oscillator, by incorporating a narrow electrical bandpass filter (Yao & 
Maleki 1996). More recently, there has been renewed interest in using this archi- 
tecture as a means for generating high-dimensional chaotic waveforms (Kouomou 
et al. 2005). 
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Figure 1. Experimental setup and corresponding mathematical block diagram of chaotic 

optoelectronic oscillator. 



The electrooptic modulator is a commercially-available lithium-niobate Mach- 
Zehnder modulator, identical to those commonly used in optical telecommunication 
systems. The input is a continuous-wave optical signal from a distributed feedback 
laser, which is split into two separate waveguide paths and then recombined, forming 
an interferometer. A voltage applied to the modulator induces an optical phase shift 
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between two arms of the interferometer through the hnear electrooptic effect. When 
the optical signals recombine, the degree to which they interfere constructively 
depends on the applied voltage. The optical power emerging from the modulator is 
then described by (Heismann et al. 1997): 

P(t) = P,cos'(l'^ + A , (2.1) 



2 K 

where Po is the continuous- wave optical power entering the modulator, v{t) is the 
voltage applied to the modulator electrodes, is the 'half- wave voltage', or the 
voltage required to produce a relative phase shift of tt between the arms of the 
interferometer, and (j)o is an angle describing the bias point of the modulator. The 
bias point is controlled either by intentionally making one arm of the interferom- 
eter longer or by adding a DC offset to the applied voltage v{t). The modulators 
described in this work had a half- wave voltage of 14 = 5.7 V and were operated at 
a bias point of (?io = — 7r/4. 

The modulator converts the applied voltage v{t) into an optical intensity modu- 
lation P{t), through the nonlinear modulation function given in equation (2.1). We 
note that this cos^(») modulation function applies to several other optical modula- 
tor structures, including liquid crystal modulators, Pockols colls (Hopf et al. 1982), 
and acoustooptic modulators (Valcc & Delisle 1985). The same nonlinearity can 
also be achieved by transmitting an electrically tunable laser through an optical 
filter that has periodic spectral transmission, such as a single-stage bircfringcnt fil- 
ter (Goedgebuer 1998) or any other single-pass interferometric filter (Blakely et al. 
2004). 

The photorcccivcr and transimpodancc amplifier produce an output voltage 
Vout{t) that is proportional to the optical power P{t), 

Voutit) = RGPit) , (2.2) 

where R is the responsivity of the photodiode (with units of A/W) and G is the 
net transimpedance gain of the system (with units of V/A.) 

The accompanying block diagram in figure 1(b) shows an equivalent mathe- 
matical diagram of the system, including the modulator, photorcccivcr, amplifiers, 
filter, and time-delayed feedback. To simplify the analysis, the voltage applied to 
the modulator is expressed in normalized units as 

xO^ff (2^3, 

and we collect all of the remaining proportionality constants into a single dimen- 
sionless factor that describes the round-trip gain of the loop, 

In terms of these dimensionlcss variables, the feedback loop relates the filter 
input r{t) to the filter output x{t) by the following nonlinear transformation and 
time delay: 

r{t)=P cos^ [x{t - r) + </>o] . (2.5) 
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For the measurements reported here, the electrical filter is a two-pole bandpass 
filter that is characterized by the linear transfer function 



His) 



(l + srL)(l + srH) ' 



(2.6) 



where tl and tu arc the time constants describing the lowpass and high-pass filters, 
respectively. In the time domain, a linear filter can be represented by state-space 
differential equations of the form: 



= Au(i) -h Br(i) 



xify = Cu{t) + Dr{t) , 



(2.7) 
(2.8) 



where r{t) is the input to the filter, x{t) is the output, u(t) is a state vector of the 
filter system, and A, B, C and D are matrices that describe the bandpass filter. For 
the two-pole bandpass filter described by equation (2.6), u(t) is a two-dimensional 
vector and the state space matrices can be expressed as 
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C=[l 0], D = 0. (2.9) 



Combining equations (2.7), (2.8), (2.9) and (2.5), the system can be described 
by the following state-space delay differential equation 



^ = Au(i) + B/3cos2 [Cu(t - t) 



00 



(2.10) 



We note that if the bandpass filter is replaced by a simple lowpass filter, then 
equation (2.10) simplifies to a scalar delay differential equation that is equivalent 
to the classic Ikeda system, originally introduced to describe bistability in optical 
cavities (Ikeda & Matsumoto 1987). 

Table 1 lists all of the parameter values used in the experiments and simulations. 
To simplify the experimental implementation, we consider here a low-frequency 
system that operates at audio frequencies, but this system can also be scaled to 
RF or microwave frequencies (Kouomou et al. 2005, Goedgebuer et al. 2002, Cohen 
et al. 2008). In practice the round-trip gain (/3) and time delay (r) were measured 
experimentally by interrupting the feedback loop at the input to the modulator and 
measuring the round-trip small-signal gain and group delay using a vector network 
analyser. The gain was controlled by varying the optical power Pq entering the 
modulator. 

In figure 2 we show calculated and measured time traces of this system, for three 
different values of the feedback strength /3, with the time delay and filter parameters 
given in Table 1. The system exhibits periodic behaviour for small values of /3, but 
the dynamics become more complex as the feedback strength is increased. Figure 3 
plots the measured and simulated bifurcation diagrams with (3 as an adjustable pa- 
rameter, showing the evolution from periodic to chaotic dynamics. Peil et al. (2009) 
reported a detailed experimental and theoretical study of the various regimes of op- 
eration of this system. In Figure 4, we plot the calculated Kaplan- Yorke dimension 
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Table 1. Experimental Parameters of System 

(The parameters used in experiments and measurements of the nonUnear chaotic oscillator. 
Here we give representative values for Po, G, R, and Vn, but in practice the factor (3 (c.f. 
equation (2.4)) was measured directly by breaking the loop and measuring the small signal, 
round-trip AC gain.) 
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Figure 2. Typical measured and calculated time traces for the nonlinear optoelectronic 
oscillator system, for feedback strengths of /3 = 1.5, 3.0, and 4.5. 

vs. f3 and vs. r, showing the progression from simple to high-dimensional chaotic 
dynamics. The Kaplan- Yorke dimension (Kaplan & Yorke 1979) was calculated 
from the spectrum of Lyapunov exponents, which were numerically computed by 
solving a linearized version of equation (2.10) (Farmer 1982). 

3. Discrete time implementation 

The optoelectronic oscillator described in §2 was introduced using a continuous- 
time delay-differential equation, but in practice, we implemented the system using 
discrete-time digital signal processing (DSP) technology. DSP provides a flexible 
platform for programmable filtering and delay operations, and offers a number of 
advantages over conventional analogue filters, especially when high-speed perfor- 
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Figure 4. Calculated Kaplan- Yorke (Lyapunov) dimension (a) as a function of the feedback 
strength /3 for fixed feedback delay r = 230 /is and (b) as a function of the the feedback 
delay r, for fixed fi — 4.5. The remaining system parameters are given in Table 1. 



mance is not required. For example, it is easy to program two digital filters to 
have identical characteristics, whereas matching of analogue filters relies on find- 
ing identical components such as resistors, capacitors and amplifiers. DSP systems 
are especially advantageous in synchronization experiments, where mismatched pa- 
rameters between nominally identical systems can otherwise impair the synchrony 
between the two systems. 

Analogue delay lines typically use either optical fibre or coaxial cables to achieve 
a time delay of L/v, where L is the length of the transmission medium and v is 
the propagation speed. Such systems cannot be scaled to large time delays because 
the required delay lines are either impractically long or prohibitively lossy. DSP 
systems, by contrast, can produce a lossless time delay that is limited only by 
the available memory and sampling rate. The use of digital processing in nonlinear 
dynamical systems dates to 1982, when Hopf et al. used a computer and ADC/DAC 
to achieve long delay times for a similar optoelectronic oscillator. Since then, digital 
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signal processing systems have dramatically improved in performance and cost, and 
are now commonplace in consumer electronics. 

Perhaps the most compelling argument in favor of DSP is that it allows real- 
time adjustment of the the gain, delay and filter coefficients - parameters that are 
typically static in analogue filter systems. In §6, we describe an adaptive control 
scheme that takes advantage of this flexibility provided by digital processing. 
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Figure 5. (a) Experimental setup showing the use of digital signal processing hardware 
to implement bandpass filter and delay, (b) Equivalent discrete-time mathematical block 
diagram of dynamical system. 



Figure 5(a) shows how the original experimental apparatus (figure 1(a)) was 
adapted to incorporate digital signal processing. The system uses the same laser, 
electrooptic modulator and photoreceiver as its continuous-time counterpart, but 
the filtering and delay are performed using a digital signal processing board. The 
DSP board uses an analogue-to-digital converter (ADC) to sample and digitize the 
input signal r{t), forming a discrete-time input sequence, 

p[n] =p{nTs) , (3.1) 

where n is an integer and Ts denotes the sampling period. The discrete-time signal 
p[n] is stored in a memory buffer to produce the desired delay and then digitally 
filtered. The output signal x[n] is then routed through a complementary digital- 
to-analogue converter (DAC) to yield the analogue output signal x{t) that drives 
the electrooptic modulator. The complete system is therefore a hybrid discrete / 
continuous-time system that retains the advantages of optical signal transmission, 
while exploiting the flexibility of discrete-time signal processing. 

The DSP board used in these experiments contains a 225 MHz floating-point 
DSP processor, 64 MB RAM, and a 16-bit ADC/DAC. The maximum sampling 
frequency was limited by the ADC/DAC chip, which was designed for audio signals. 
Except where noted, we used a sampling rate of l/T^ = 96 kS/s in these experi- 
ments, although the system could be scaled to higher frequencies by replacing the 
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ADC/DAC hardware. The lowpass filter in the feedback loop restricts dynami- 
cal behaviour to frequencies well below the Nyquist frequency, ensuring that the 
sampling docs not contribute signficantly to the filtering. Higher-performance field- 
programmable gate array (FPGA) boards could perform the same operations at 
sampling rates as high as 1 GS/s. 

Figure 5(b) shows a mathematical block diagram of the discrete-time system. 
In this system, k denotes the feedback delay, which we take to be an integer num- 
ber of timesteps, and the dynamical filter is described by discrete difference equa- 
tions rather than differential equations. The digital filter was designed to act as 
a two-pole bandpass filter that approximates the response of the continuous-time 
filter described in equations (2.6)-(2.8). The discrete-time transfer function H{z) 
is obtained from the continuous-time transfer function H(s) by applying a bilinear 
transform with frequency pre- warping (Oppenheim et al. 1999). This process yields 
the following equivalent discrete-time transfer function 



1 



(1 







if(.) = -(!-. .)(1 + .h)-_^^^_,^^^_^^^_,^ 



(3.2) 



where zl and zh are the poles of the discrete-time filter, which are related to the 
time constants tl and th and sampling period Tg by 



Zh = 



1 — tan 



ZL = 



2tl 



1 -|- tan 



T 

2^ 



(3.3) 



The discrete-time filter can be represented by state-space evolution equations 
analogous to equations (2.7) and (2.8), 



u[n -I- 1] = Au[n] + Br[n] 
x[n] = Cu[n] -|- Dr[n] 



(3.4) 
(3.5) 



where r[n] is the filter input, x[n] is the output, and u[n] is a two-dimensional state 
vector. For the filter described in equation (3.2), the state space matrices can be 
expressed as 
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{zL + Zh) -Zl 
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{1 - zl){1 + zh){1 + zlZh) 
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1 



D = -{1-zl){1 + zh). (3.7) 



AzlZh 

The filter input is related to the filter output through a nonlinearity and delay, 

r[n] = fi cos^ {x[n — k] + 0o) > (3. 



where fi is the round- trip gain defined in equation (2.4) and the delay is chosen to 
be fc = 22, which, at a sampling rate of 96 kS/s, corresponds to a feedback delay 
of 830 jjis. 
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4. Coupled systems and synchronization 

An interesting property of chaotic systems is that two systems, when properly 
coupled together, can synchronize with one another and evolve along the same 
chaotic orbit (Fujisaka & Yamada 1983, Pecora & Carroll 1990, Pikovsky et al. 
2001, Boccaletti 2008). Many proposed applications of chaos, including secure com- 
munication systems, sensor networks, and data assimilation and prediction, rely 
on this phenomenon of synchronization between chaotic oscillators (Kanter et al. 
2008, Argyris et al. 2005, Golubitsky et al. 2005, Boccaletti et al. 2006, Arenas et 
al. 2008). There have been some analytical studies of the coupling threshold re- 
quired for synchronization in delayed-feedback systems (Pyragas 1998, Biinner & 
Just 1998). Peil et al. (2007) reported some experimental measurements and the- 
oretical models of synchronization between time-delayed optoelectronic oscillators 
like those discussed here. We seek in the this section to more thoroughly investigate 
how two such systems can be coupled together, and the conditions under which they 
can synchronize. 
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Figure 6. (a) Block diagram of two linearly-coupled optoelectronic chaotic oscillators, 
where the coupling and delays are taken to be in the optical path connecting the two 
systems, (b) Equivalent system, obtained by commuting the coupling and delay with 
the bandpass filter. In practice, the coupling is implemented optically as in (a), but for 
convenience, we analyse the equivalent scenario depicted in (b). 



The block diagram in figure 6(a) shows the most general type of linear optical 
coupling between two systems. In this case, we imagine that the optical signal 
emerging from the modulator in system 1 is split and fed back into both systems. 
The constants Pn and rn denote the self-feedback gain and delay for system 1, and 
/3i2 and T12 describe the coupling from system 1 — > 2. Similarly, P22, and r22 are the 
self-feedback parameters of system 2, and P21 and T21 describe the coupling from 
2^1. We assume that the bandpass filters between the two systems are identical. 

The filter (7J(s)), gain iPij) and delay (r^j) are all linear, time-invariant oper- 
ations, and they can therefore be freely permuted without changing the dynamics 
of the system. Using these arguements, one can transform the optically coupled 
system shown in figure 6a to the equivalent system shown in figure 6b, where the 
coupling instead applies to the electrical signals Xj it) emerging from the bandpass 
filters. This coupling configuration can be described by the following coupled delay 
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differential equations: 

^ = Aui(i) + B cos^ [/3iiCui(i - ni) + /32iCu2(t - T21) + cPi] (4.1) 

^ = Au2(i) + B cos^ [/322Cu2(i - T22) + /3i2Cui(t - ns) + , (4.2) 

where ui and U2 are the state-vectors for the bandpass filters in oscillators 1 and 

2, respectively. 

To understand the conditions under which synchrony can occur, we begin by 
assuming that a synchronous solution exists, 

Ui{t)=U2it-To) = u{t), (4.3) 

where we have allowed for lag synchrony with a time delay tq. Upon substituting 

this assumption into equations (4.1) and (4.2), wc obtain self-consistent dynamical 
equations for u(t) only under the following conditions: 

<t>l=4>2, /3ll = fe, /3l2=/321, Tii=r22, To = -(t2i - T12) . (4.4) 

While these conditions arc necessary for a synchronous solution to exist, they do 
not guarantee the stability of this solution. 

We now further restrict our attention to cases in which the systems synchronize 
in a state that obeys the same dynamical equation as that of an uncoupled, isolated 
system described by parameters /3, r and 4>q. This lifts the constraint that = /322 
and /3i2 = , but imposes the following additional conditions for synchrony: 

011 + 021 = 022 +012=0 (4.5) 
Til = T2I = T22 =Ti2=T (4.6) 

(l)i=4>2= (t>o (4.7) 
To = . (4.8) 

In this scenario, which is termed 'diffusive coupling', the constraint on the coupling 
conditions (equation (4.5)) can be cast in terms of two dimensionless parameters 
Ki and K2-, defined through the relations 

021 = 011 = (1 - «;i)/3 (4.9) 

012 = K20, /322 = (1 - «2)/3 . (4.10) 

With this definition, (1 — Ki) and Ki describe relative proportions of self- feedback 
vs. cross-coupled feedback, respectively, entering system 1, and K2 has a similar 
interpretation for system 2. 

Figure 7 presents the block diagram of two diffusively- coupled oscillators. In 
order to make the equations comparable in form to the single-oscillator system 
described in §2, wc have factored out a common scale factor from all four of the 
coupling terms and commuted this scale factor with the bandpass filter H[s). 

The diffusively coupled oscillator system shown in figure 7 is described by the 
following coupled equations 

^ = Aui(t)+B/3cos2(c[(l-Ki)ui(t-T) + MU2(t-r)] +<^o) (4.11) 
^ = Au2(0+B/3cos2(c[(l-K2)u2(t-r)+K2Ui(t-r)] + ./.q) . (4.12) 
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Figure 7. Block diagram of two diffusively-coupled oscillators. The coupling is constructed 
in a way that guarantees that the resulting system admits a synchronous solution of the 
form xi{t) = X2{t) = x{t) where x{t) exhibits the same dynamical behaviour as that of an 
isolated system. 



These equations can be seen to admit an isochronally synchronized solution that, 
when synchronized, satisfies the same equation (2.10) given earlier for an isolated 
system. 

To investigate the stability of the synchronized solution, wc perform following 
change of variables 

u+W-^[uiW + U2W], U-W-^[uiW-U2W] , (4.13) 

where the difference u_ (<) is expected to converge to zero for a stable synchronous 
solution. Expressing equations (4.11) and (4.12) in terms of the u±, and linearizing 
about the synchronous state, we find 



= Au+(t) + B/3cos^ Cu+(t - T 



(4.14) 



Au_(i) +B/3sin(2Cu+(t - t) + 20o)(ki + K2 - l)u-(i - r) . (4.15) 



du^ 

"dt 
du_ 

"dt 

Comparing equations (4.14) and (2.10), we see that VL+{t) satisfies the same 
dynamical equation as isolated system, as expected. The two coupling parameters 
appear in the second equation only in the combination [ki + K2). We therefore 
conclude that for a given /3, r and (pQ, the stability of the synchronous solution 
depends only on the sum (ki + ^2), but not on the values of ki and K2 individually. 

Furthermore, in the special case that ki + K2 = 1, equation (4.15) simplifies to 



(iu_ 
dt 



Au_(t) 



(4.16) 



Because the linear bandpass filter is stable (i.e., A has negative eigenvalues) the 
difference vector u_(t) will always converge to zero according to the filter time- 
constants tl and th whenever ki -I- K2 = 1. 

Figure 8 plots the measured and simulated normalized root-mean square (RMS) 
synchronization error as a function of the coupling strength k for the case of bidi- 
rectional symmetric coupling, i.e., ki = K2 = showing the regimes in which the 
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Figure 8. Measured and simulated normalized synchronization error as a function of k for 
/3 = 6 and /? = 8. The measurements and simulations were conducted using symmetric 
bidirectional coupling, (ki = K2 = h). 

two systems synchronize. We define the normalized synchronization error as 



where (•) indicates a time average. The normalized error Gx is zero in the case 
of a synchronized solution, but approaches 1 in the limit that the two signals are 
identically distributed, but uncorrelated. While the experimental measurements 
and simulations were performed by taking ki = K2, the results can be generalized 
to other combinations of ni and k,2 because the synchronization condition depends 
only on the sum ki + K2. 

As indicated in figure 8, the two systems synchronize unconditionally for the 
special case that ki + K2 — 1, and they synchronize for a range of (ki + K2) cen- 
tered symmetrically about this point. The range of values over which the system 
synchronize is found to depend on the feedback gain f3. In general, we observed 
that the higher values of /3 (and higher Lyapunov dimension) yield a narrower 
synchronization regime. 



In addition to knowing whether two systems synchronize, it is also important to un- 
derstand the rate at which they converge to a synchronous state, which is quantified 
by the transverse Lyapunov exponent (Fujisaka & Yamada 1983, Pecora & Carroll 
1998). The transverse Lyapunov exponent (TLE), denoted A^, defines the average 
exponential rate at which a pair of coupled identical oscillators converge or diverge 
in phase space. A negative TLE corresponds to converging trajectories, indicating 

Article submitted to Royal Society 




(4.17) 



5. Synchronization — Transient dynamics 



14 



T. E. Murphy and others 



stable synchronization, while a positive exponent indicates diverging trajectories 
that do not synchronize. 

The TLE defines an important timcscalc in applications such as chaotic com- 
munication and sensor networks that rely on synchronization. In chaotic sensor 
networks with time-varying coupling, the TLE limits the speed of perturbations 
that the system can track. In a chaotic communication system, the TLE limits the 
attainable bit-rate that can be successfully decoded. 

Equally important is the (positive) maximal Lyapunov exponent, which de- 
scribes the rate at which initially synchronous solutions diverge from one another 
when they are decoupled. This divergence rate is important in data assimilation and 
prediction applications, which use synchronization to predict the future behaviour 
of a dynamical system. Here we present a numerical and experimental study of the 
transient synchronization and desynchronization dynamics of two coupled chaotic 
optoelectronic oscillators. 

One method to determine the transverse Lyapunov exponent is to suddenly cou- 
ple two independent and identical chaotic oscillators. By analysing the transition 
from the initially uncorrelated dynamics to a synchronous state, we can determine 
the (finite time) transverse Lyapunov exponent of the system. Conversely, if the 
two systems are initially synchronized, the coupling can be suddenly turned off, 
allowing the trajectories to exponentially diverge. By measuring the rate of expo- 
nential divergence, we find the maximal Lyapunov exponent of the system (Cohen 
et al. 2008). Unlike conventional methods, which require numerical solution of a 
linearized system of equations, this approach can be applied even in cases when an 
exact model of the physical system is unavailable or impractical. As long as two 
experimental systems can be made to synchronize, the Lyapunov exponents de- 
scribing synchronization and desynchronization can be determined from transient 
time-series analysis. 

This method of determining the Lyapunov exponent is illustrated in Figure 9, 
which shows the exponential convergence and divergence of two coupled chaotic 
optoelectronic oscillators. In figure 9(a)-(b), the two oscillators were initially un- 
coupled for i < 0, but the coupling was suddenly enabled at t = 0. Specifically, 
for t > the systems were bidirectionally coupled as shown in figure 7 with 
Ki = K2 = K = 0.4375. Figure 9(a) plots the measured outputs xi{t) and X2{t) 
for one representative case, showing the transition from uncorrelated to synchro- 
nized dynamics. Figure 9(b) shows the absolute difference \xi{t) — x^it)], smoothed 
with a 100 fis sliding window average, and plotted on semilogarithmic axes to clearly 
show the exponential convergence. By fitting an exponential relation to this curve, 
we determine the (negative) transverse Lyapunov exponent. Figure 9(c)-(d) show 
similar data obtained when two initially synchronized systems are decoupled at 
t = 0, allowing them to exponentially diverge. In this case, the (positive) maximum 
Lyapunov exponent Ai is similarly determined by finding the best-fitted slope to 
the smoothed logarithmic difference between the two traces. 

When determining the Lyapunov exponent using this method, the exponential 
convergence or divergence is estimated only over a finite fitting interval T. In prac- 
tice, the allowable fitting interval is restricted by the synchronization error floor, 
which is caused by noise and mismatches between the two systems (Shahverdiev 
et al. 2005). In numerical simulations, the convergence/ divergence can be observed 
over many orders of magnitude, and we can therefore fit the exponential relation 
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Figure 9. (a) Experimentally measured time series showing synchronization of two coupled 
chaotic oscillators. The two systems were uncoupled for t < and symmetric bidirection 
coupling was abruptly enabled at t = 0. (b) Measured absolute difference \x\{t) —- X2{t)\ 
plotted on a logarithmic scale, and smoothed to show exponential convergence of trajec- 
tories. By fitting a line to this slope, one can estimate the finite-time transverse Lyapunov 
exponent At, which characterizes the timescale over which synchronization occurs, (c) 
Experimentally measured time series showing divergence of two initially synchronized sys- 
tems, when the coupling is disabled at t = 0. (d) The finite-time maximal Lyapunov 
exponent Ai is estimated by measuring the average exponential divergence rate. 



over a larger time window T . If the fitting window is long enough to span the en- 
tire chaotic attractor, this calculation reveals the 'global' or 'asymptotic' Lyapunov 
exponent. For a short fitting interval, the trajectory remains only in a localized 
portion of the chaotic attractor, and thus we obtain only a 'local' Lyaponov expo- 
nent. The local Lyapunov exponents vary about an attractor, and their statistical 
distribution depends upon the dynamical nature of the coupled system. 

In figure 10, we show distributions of local Lyapunov exponents for three choices 
of fitting time T for 10^ simulated time-series. The histograms labeled (a) show the 
distribution of transverse local Lyapunov exponents, obtained by simulating two 
initially independent systems that are suddenly coupled together with ki = K2 = 0.4 
'At t — Q. The histograms labeled (b) histograms show the distribution of maximum 
Lyapunov exponents, obtained by simulating two initially synchronized systems 
that are suddenly released at t = 0. In all cases, the histogram is Gaussian near its 
peak and has non-Gaussian tails. The mean of each distribution converges to the 
global or average (transverse) Lyapunov exponent A, and the standard deviation 
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Figure 10. Histogram showing the distribution of finite-time transverse Lyapunov ex- 
ponents, measured over time intervals of 2, 4, and 8 ms. The transverse and maximal 
Lyapunov exponents were determined by numerically simulating the coupled system and 
fitting the convergence or diverence to an exponential relation, as depicted in figure 9(b) 
and (d). 

narrows in proportion to T"^/^, as expected (Prasad & Ramaswamy 1999, Ott 
1993). 

When a mathematical model of the system is available, the transverse Lyapunov 
exponents can also be calculated using the master stability function technique (Fu- 
jisaka & Yamada 1983, Pecora & Carroll 1998); i.e., by linearization about the 
synchronized chaotic solution. Figure 11 compares the distribution of local trans- 
verse Lyapunov exponents obtained using both methods. In figure 11(a) we plot (in 
grayscale) the distribution of local transverse Lyapunov exponents as a function 
of the coupling strength n {— ki = K2). These histograms were obtained using 
time-series analysis to estimate the exponential convergence, as illustrated in figure 
9(b). Because the time traces were initially uncorrelated, this method only applies 
when the TLE is negative, corresponding to convergent time series. As anticipated 
from equation (4.15), the systems converge unconditionally when ki = K2 = 0.5. 
Figure 11(b) plots the same distribution of TLEs, obtained by numerically solving 
the linearized system of equations. Here, the linearized equations are sensitive to 
both positive and negative phase space growth, so the distributions can go above 
zero. Apart from this expected difference, the correspondence between these two 
methods is remarkably good. 

In figure 11(c), we plot the transverse Lyapunov exponents obtained from ex- 
perimentally measured converging time-series of coupled systems. At each value of 
AC, we measured the convergence rates At for 100 pairs of time-series. The mean 
At and standard deviation obtained by fitting the data to a Gaussian distribution 
are shown as the dots and bars on the figure respectively. For comparison, the line 
indicates the 'global' TLE computed using the linearized system of equations. The 
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Figure 11. Distribution of transverse Lyapunov exponents as a function of the coupling 
strength k, determined (a) by time-series analysis of converging transients and (b) by 
directly solving the linearized transverse equations of the coupled system, (c) Comparison 
of calculated global transverse Lyapunov exponent and the measured finite-time transverse 
Lyapunov exponent, determined by measuring the transient convergence. The data points 
and error bars indicate the average and standard deviation of the measured statistical 
distribution of At, using a finite fitting time of 4 ms. 

experimental data agree well with the numerical simulations, which demonstrates 
that time-series analysis of converging experimental signals is a powerful technique 
for quantifying the transverse Lyapunov exponents of a system, even if a numerical 
model were unavailable. 

6. Adaptive synchronization 

As shown in §5, synchronization can depend on the coupling between oscillators. In 
a practical network consisting of spatially separated chaotic oscillators, time- varying 
environmental conditions can cause the coupling to vary unpredictably. Then, in 
order to maintain synchronization it is essential to dynamically compensate for these 
variations. Recently, several algorithms have been developed to maintain or produce 
synchrony in a network of chaotic oscillators (Zhou & Kurths 2006, Ito & Kaneko 
2001, De LeUis 2008, Feki 2003). Sorrentino & Ott (2008, 2009) proposed and 
simulated an adaptive algorithm to estimate and track a priori unknown coupling 
changes in a network of chaotic oscillators. The estimate is then used to compensate 
for the environmental perturbations thereby ensuring synchrony. In this section, we 
present an experimental demonstration of this scheme using a pair of nonlinear time- 
delayed optoelectronic feedback loops described in §2. A DSP board, incorporated as 
part of the feedback loop (see figure 5), enables us to perform real time computations 
allowing the implementation of the adaptive tracking algorithm. 

In our experimental setup, shown in figure 12(a), we consider two optoelectronic 
feedback loops that are unidirectionally coupled through a time- varying communi- 
cation channel, which is described by a coupling factor K{t). An adaptive control 
scheme is implemented in the receiver in order to maintain synchrony between the 
two systems and, in the process, determine an estimate of the channel condition. 

Figure 12(b) shows an equivalent discrete-time mathematical block diagram of 
the two unidirectionally coupled systems with a time-varying channel. Here we de- 
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Figure 12. (a) Experimental setup of unidirectionally coupled feedback loops, in which 
the coupling factor k is allowed to vary slowly, to simulate the effect of an atmospheric 
perturbation or environmental disturbance, (b) Equivalent discrete-time block diagram of 
two oscillators, unidirectionally coupled over a time-varying channel. The receiver has no 
prior knowledge of «:[n], and must therefore form an estimate, denoted in order to 
keep the two systems sychronized. 



note the channel coupling by K[n]. The receiving system has no a priori knowledge 
of and must therefore form an estimate, denoted K[n], in order to maintain 
isochronal synchrony. As before, the discrete-time bandpass filters H{z) are gov- 
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erned by the state-space equations 

u^[n+^] = Aui[n]+Bri[n] (6.1) 
Xi[n] = Cu,[n] + Dri[n\ (6.2) 
(i = l,2), 

where ri[n], {i = 1,2) are the filter inputs, Xi[n] are the corresponding filter outputs, 
and A, B, C and D describe the bandpass filter. The filter outputs are fed back to 
the inputs through a nonlinearity and time delay according to 

n [n] = P cos^(a;i [n-k] + 0o) (6.3) 
r2[n] = /3cos^[(l — /t[n — A;])a;2[n — k] + K[n — k]xi[n — k]+ (f>o] , (6.4) 

where K[n] is the local estimate of the channel coupling. 

One can clearly sec that these equations admit a synchronous solution in the case 
that K[n] = K[n], i.e., provided the receiver tracks the coupling strength K[n\. The 
analysis presented in §4 showed that for static coupling, i.e., when K[n] = k[ti] = k, 
the synchronous solution is stable over a continuous range of values of k. This result 
suggests that if K[n] varies slowly, while remaining within the bounds required for 
synchrony stability, the systems could stay synchronized as long as the receiver is 
able to track the variation with sufficient accuracy. We emphasize that the receiver 
does not have direct knowledge of K[n\ but only receives the product K[n]a:i[n], as 
shown in figure 12(b). 

Sorrentino & Ott (2008, 2009) prescribed a strategy in which the local factor 
K[n] is adjusted in a way that minimizes the average synchronization error. This 
yields the following estimate /t[n], 



where (•)lpf denotes an exponentially- weighted moving average, wliic;li is equiva- 
lent to a discrete-time low-pass filter. This averaging process can be implemented 
with the following discrete-time iterative equations: 

N[n] = ZoN[n - 1] + (1 - zo)K[n]xi[n]x2[n] (6.6) 
D[n] = zoD[n - 1] + (1 - zo)xl[n] , (6.7) 

where the forgetting factor Zq is the pole of the discrete-time low-pass filter. The 
time- window over which the averaging is performed is approximately Ts{l — Zo)~^ , 
where Ts is the sampling period. We note that, as required, the adaptive scheme 
described by equation (6.5) relies only on the product K;[n]2;i[n] and X2[n] to form 
the estimate K[n]. In a high-speed application, the lowpass filter could easily be 
implemented using an electrical mixer in place of discrete-time averaging filter. 

We experimentally demonstrated the adaptive synchronization scheme using 
a pair of coupled nonlinear optoelectronic oscillators, as shown in figure 12(a). 
For these experiments, the bandpass filters were adjusted to have a passband of 
100 Hz - 2.5 kHz, the DSP sampling frequency was reduced to 24 kS/s, and the 
time delay was measured to be fc = 36 timesteps, or 1.5 ms. We chose a feedback 
strength of /3 = 3.58, which, under these conditions, was found to yield robust 
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chaotic behaviour. The lowpass fiher used in the adaptive synchronization rule was 
implemented with a forgetting factor of zq = 0.95, which corresponds to a filter 
response time of 208 /iS. 
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Figure 13. (a) Measured and (b) simulated response of adaptive coupling system to a 
sudden change in In these experiments, the coupling strength n was changed abruptly 
from 0.80 to 1.13 at t = 0. The adaptive synchronization scheme automatically adjusts R 
in response to this variation. Here we plot both the tracking signal R{t) and the difference 
xi{t) — X2{t), showing the initial loss of synchrony followed by recovery. 



Figure 13 presents experimental measurements and numerical simulations show- 
ing how both the synchronization error and tracking signal R[n\ respond to an 
abrupt change in the coupling from k — 0.8 to k = 1.13. For t < Q the coupling 
strength k was held constant at k = 0.80. Under these conditions the receiver forms 
the correct estimate R = 0.8, which gives a small synchronization error. At < = the 
coupling strength was switched abruptly to k = 1.13, which causes the two loops 
to briefly lose synchrony. However, the receiver adaptively readjusts the parameter 
K to track K[n] and the synchrony is regained. The numerical simulations shown in 
figure 13(c) and (d) exhibit similar behaviour. The response time of the adaptive 
synchronization method was found to be limited primarily by the exponentially- 
weighted moving average filter. In a separate work we studied the ability of this 
adaptive scheme to track sinusoidal variations in coupling, and we quantified the 
limitations on the magnitude and frequency of the perturbation that can be tracked 
(Ravoori et al. 2009). 
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7. Conclusion 

Many potential applications such as secure communication, sensor networks, spread- 
spectrum communication, chaotic radars and random number generators could ben- 
efit from a nonlinear dynamical system that is simple to model, easy to implement, 
and capable of generating robust, high-dimensional, chaotic waveforms. This paper 
presents a comprehensive analysis and characterization of a nonlinear optoelectronic 
feedback system that meets these criteria. The system uses elcctrooptic modulation 
and optical transmission, and it can therefore take advantage of the vast array of 
low-cost, high-speed, widely available components originally developed for fibre- 
optic communication networks. Wc describe a new approach in which the delayed 
electrical feedback and filtering is implemented using real-time digital signal pro- 
cessing. This greatly facilitates matching of filter characteristics between systems, 
and also allows for real-time control and adjustment of the feedback parameters 
something that could not be easily accomplished with traditional analogue signal 
processing. 

Because most of the aforementioned applications of chaotic signals require syn- 
chronization between two or more systems, we explore the conditions under which 
coupled system will synchronize. We present a new technique to experimentally 
quantify the rate of convergence when two systems are coupled and the rate of 
divergence when they are released. Finally, we demonstrate an adaptive technique 
that automatically maintains synchronization between coupled systems, in the pres- 
ence of an unknown and time- varying coupling between the two. 

This work was supported by DOD MURI grant (ONR N000140710734) and the US-Israel 
Binational Science Foundation. 
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